SUBROUTINE vfunapprox_thetav_q(emaxin,betain)
  USE Commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'

  REAL(8), INTENT(IN)  :: emaxin(Ntheta_vot,Nwealth,Nq)
  REAL(8), INTENT(OUT) :: betain(MMtheta_vot_steal)
  REAL(8), ALLOCATABLE :: y(:), xvar(:,:), res(:,:,:), vfunout(:,:,:), ones(:), betastart(:), xpx(:,:), xpy(:), invxx(:,:), res_final(:), work(:)
  REAL(8)              :: Rsquared, av_res(Nwealth), av_res_temp(Nwealth), v_slope, v_slope_p
  INTEGER              :: ii, i, j, iwealth, itheta, info, MM_t, lwork, iq
  CHARACTER :: uplo='U'

  ii = Ntheta_vot*Nwealth*Nq

  ALLOCATE(y(ii), res(Ntheta_vot,Nwealth,Nq), ones(ii), vfunout(Ntheta_vot,Nwealth,Nq), res_final(ii))
  y = 0d0; res = 0d0

  MM_t = 1d0 + Npol_theta_vot + Npol_q
  ALLOCATE(xvar(ii,MM_t),xpx(MM_t,MM_t), xpy(MM_t), invxx(MM_t,MM_t), betastart(MM_t))

  xpx = 0d0; xpy = 0d0
  betastart = 0d0
  betain = 0d0
  !Average the residuals for each iwealth
  av_res_temp = 0d0
  DO iwealth = 1, Nwealth
     DO itheta = 1, Ntheta_vot
        DO iq = 1, Nq
           av_res_temp(iwealth) = av_res_temp(iwealth) + emaxin(itheta,iwealth,iq)
        END DO
     END DO
  END DO
  av_res_temp = av_res_temp/DBLE(Ntheta_vot*Nq)

!!$  IF (flag_here == 2858) THEN
!!$     WRITE(*,'(a,10f15.7)') 'av_res_temp',  av_res_temp(:)
!!$     DO iq = 1, Nq
!!$        WRITE(*,'(a,1i4,10f15.7)') 'emaxin', iq, emaxin(iq,:)
!!$     END DO
!!$  END IF

  !Compute the residuals V - (alpha_1+alpha_2*q+alpha_3*q^2)
  DO iwealth = 1, Nwealth
     DO itheta = 1, Ntheta_vot
        DO iq = 1, Nq

           res(itheta,iwealth,iq) = emaxin(itheta,iwealth,iq) - av_res_temp(iwealth)
!!$        IF (myrank == 0) WRITE(*,'(a,2i4,3f15.7)') 'res', itheta, iwealth, res(itheta,iwealth), emaxin(itheta,iwealth), av_res_temp(iwealth)

        END DO
     END DO
  END DO

  !We compute the coefficients on public consumption in V = f(sav)+alpha_1+alpha_2*q+alpha_3*q^2
  ones = 1d0
  i = 0
  DO iwealth = 1, Nwealth
     DO itheta = 1, Ntheta_vot
        DO iq = 1, Nq

           i = i + 1
           y(i) = res(itheta,iwealth,iq)
           xvar(i,1:MM_t) = (/ ones(i), theta(itheta), qconsvec_g(iq), DBLE(itheta-1)*qconsvec_g(iq) /)

!!$             IF (flag_here == 260) write(*,'(a,3i4,3f14.8)') 'y inside', itheta,iwealth,iq,y(i),emaxin(itheta,iwealth,iq),av_res_temp(iwealth)

        END DO
     END DO
  END DO
  
  xpx = MATMUL(TRANSPOSE(xvar(:,:)),xvar(:,:))/SIZE(y)
  xpy = MATMUL(TRANSPOSE(xvar(:,:)),y(:))/SIZE(y)
     
  invxx = 0d0
  DO j=1, MM_t
     invxx(j,j)=1d0
  END DO

  LWORK = MM_t
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))  
  CALL DPOSV(uplo,MM_t,MM_t,xpx,MM_t,invxx,MM_t,info)
  IF (info /= 0) THEN
     DO j = 1, ii
        WRITE(*,'(a,1i3,7f12.6)') 'y,x',j,y(j),xvar(j,1:MM_t)
     ENDDO
     WRITE(*,*) 'MATRIX COULD NOT INVERT vfun_thetav ', 'info', info
  END IF
  
  betastart(:) = MATMUL(invxx,xpy)

  IF (flag_Tend == 1) THEN
     v_slope_p = (av_res_temp(2)-av_res_temp(1))/(wealth(2) - wealth(1))
     av_res(Nwealth) = av_res_temp(Nwealth)
     DO iwealth = Nwealth-1, 1, -1
        IF (av_res_temp(iwealth+1) - av_res_temp(iwealth) > small_no) THEN
           av_res(iwealth) = av_res_temp(iwealth)
           V_slope = (av_res_temp(iwealth+1)-av_res_temp(iwealth))/(wealth(iwealth+1) - wealth(iwealth))
        ELSE        
           V_slope = inc_slope*V_slope_p+V_slope_p
           av_res(iwealth) = V_slope*(wealth(iwealth) - wealth(iwealth+1)) + av_res(iwealth+1)
        END IF
        V_slope_p = V_slope
     END DO
  ELSE
     av_res = av_res_temp
  END IF
!!$  betain(:) = (/ av_res(:), betastart(:) /)
  betain(1:Nwealth+MM_t-Npol_q) = (/ av_res(:), betastart(1:MM_t-Npol_q) /)
  betain(Nwealth+MM_t-Npol_q+Npol_st+1:Nwealth+MM_t+Npol_st) = betastart(MM_t-Npol_q+1:MM_t)

!!$  IF (flag_here == 260) THEN
!!$     write(*,'(a,18f14.8)') 'betain try', betain
!!$     write(*,'(a,18f14.8)') 'betastart', betastart
!!$  END IF
  
  !Compare the actual value function with the approximation
  i = 0
  DO iwealth = 1, Nwealth
     DO itheta = 1, Ntheta_vot
        DO iq = 1, Nq

           i = i + 1
           vfunout(itheta,iwealth,iq) = av_res(iwealth) + betastart(1) + betastart(2)*theta(itheta) + betastart(3)*qconsvec_g(iq) + betastart(4)*DBLE(itheta-1)*qconsvec_g(iq)
           res_final(i) = emaxin(itheta,iwealth,iq) - vfunout(itheta,iwealth,iq)
           y(i) = emaxin(itheta,iwealth,iq)

        END DO
     END DO
  END DO

  Rsquared = (SUM((y-SUM(y)/SIZE(y))**2d0) - SUM(res_final**2d0)) / SUM((y-SUM(y)/SIZE(y))**2d0)

  IF (Rsquared < 0.75d0 .AND. flag_Tend == 0) THEN
!!$  IF (Rsquared < 0.75d0) THEN
!!$  IF (flag_here == 1) THEN
!!$  IF (flag_Tend == 1) THEN
     IF (myrank == 0) THEN
        write(*,*) 'Rsquared vfunapprox_thetav_q', Rsquared
        DO iwealth = 1, Nwealth
           DO itheta = 1, Ntheta_vot
              DO iq = 1, Nq

                 write(*,'(a,3i4,5f20.10)') 'vfunapprox_thetav_q', itheta,iwealth,iq,emaxin(itheta,iwealth,iq), vfunout(itheta,iwealth,iq), av_res(iwealth) + betastart(1) + betastart(2)*theta(itheta) + betastart(3)*qconsvec_g(iq) + betastart(4)*DBLE(itheta-1)*qconsvec_g(iq), av_res_temp(iwealth), betastart(1) + betastart(2)*theta(itheta) + betastart(3)*qconsvec_g(iq) + betastart(4)*DBLE(itheta-1)*qconsvec_g(iq)

              END DO
           END DO
        END DO
        write(*,*) 'betas', betain(:)
     END IF
     flag_here = 146
  END IF
  
  DEALLOCATE(y,xvar,res,vfunout)

END SUBROUTINE vfunapprox_thetav_q
